Chapter 4 Community composition

load("data/data.Rdata")

4.1 Taxonomy overview

4.1.1 Stacked barplot

# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  filter(count > 0) #filter 0 counts

# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)

ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ type,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
  facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label

4.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_qdedwpfjwj83upocimp0
phylum mean sd
p__Bacteroidota 0.376352935 0.205078492
p__Bacillota_A 0.250845368 0.165219604
p__Bacillota 0.116698680 0.146369286
p__Pseudomonadota 0.095767942 0.159976047
p__Campylobacterota 0.052429876 0.092626270
p__Verrucomicrobiota 0.029518627 0.072768229
p__Desulfobacterota 0.023304680 0.036356268
p__Actinomycetota 0.012917195 0.107880476
p__Chlamydiota 0.010449046 0.059349062
p__Fusobacteriota 0.010359534 0.028167393
p__Cyanobacteriota 0.009098787 0.016424861
p__Bacillota_C 0.004658039 0.006626014
p__Spirochaetota 0.003962784 0.012243051
p__Bacillota_B 0.002436629 0.004905803
p__Elusimicrobiota 0.001199879 0.006464719 
phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal() + 
        theme(legend.position="none") +
        labs(y="Phylum",x="Relative abundance")

4.2 Taxonomy boxplot

4.2.1 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_g0dnrymakwh9vc5ughs4
family mean sd
f__Bacteroidaceae 2.189238e-01 0.1404298747
f__Lachnospiraceae 1.389745e-01 0.1089608538
f__Tannerellaceae 1.016053e-01 0.0802366634
f__Helicobacteraceae 5.197482e-02 0.0922005028
f__Mycoplasmoidaceae 3.651280e-02 0.0753015455
f__Erysipelotrichaceae 3.461428e-02 0.0450141949
f__UBA3700 3.378612e-02 0.0555421417
f__Marinifilaceae 2.741603e-02 0.0272043836
f__ 2.740335e-02 0.0833966811
f__DTU072 2.668734e-02 0.0526713763
f__Rikenellaceae 2.665357e-02 0.0462652505
f__Enterobacteriaceae 2.660196e-02 0.0913429878
f__Coprobacillaceae 2.521928e-02 0.0887342893
f__Desulfovibrionaceae 2.330468e-02 0.0363562684
f__Ruminococcaceae 1.838184e-02 0.0427420605
f__LL51 1.738859e-02 0.0683714736
f__Rhizobiaceae 1.578130e-02 0.0766280394
f__UBA3830 1.560351e-02 0.0451659123
f__Mycobacteriaceae 1.228390e-02 0.1079320546
f__Akkermansiaceae 1.213003e-02 0.0315247341
f__Chlamydiaceae 1.044905e-02 0.0593490616
f__Fusobacteriaceae 1.035953e-02 0.0281673926
f__CAG-239 9.034415e-03 0.0150424841
f__Enterococcaceae 8.322037e-03 0.0463917038
f__Gastranaerophilaceae 7.638375e-03 0.0143794214
f__Oscillospiraceae 6.565519e-03 0.0077974910
f__UBA1997 6.303806e-03 0.0307918484
f__Streptococcaceae 6.291980e-03 0.0340225424
f__UBA1242 4.618700e-03 0.0152906152
f__Brevinemataceae 3.962784e-03 0.0122430505
f__Acutalibacteraceae 3.335081e-03 0.0108974502
f__UBA660 3.159125e-03 0.0138778160
f__Clostridiaceae 2.808963e-03 0.0170355922
f__RUG11792 2.784774e-03 0.0249664653
f__Peptococcaceae 2.436629e-03 0.0049058034
f__MGBC116941 2.134904e-03 0.0093302936
f__Acidaminococcaceae 1.921070e-03 0.0050068934
f__CAG-508 1.786832e-03 0.0063826540
f__Anaerovoracaceae 1.551174e-03 0.0036295746
f__Moraxellaceae 1.468042e-03 0.0096884882
f__RUG14156 1.460412e-03 0.0045605184
f__Staphylococcaceae 1.345783e-03 0.0050584305
f__Elusimicrobiaceae 1.199879e-03 0.0064647185
f__CAG-288 9.379861e-04 0.0059852367
f__Anaerotignaceae 8.884602e-04 0.0040242420
f__CALVMC01 7.428930e-04 0.0043360033
f__Eggerthellaceae 6.332937e-04 0.0021152288
f__Massilibacillaceae 6.162148e-04 0.0016267915
f__UBA1820 4.680638e-04 0.0012990360
f__Arcobacteraceae 4.550604e-04 0.0050029784
f__CAG-274 4.466884e-04 0.0021903937
f__Burkholderiaceae_C 3.656162e-04 0.0047810524
f__Muribaculaceae 3.575580e-04 0.0009726557
f__UBA932 3.379555e-04 0.0011520742
f__Hepatoplasmataceae 2.954146e-04 0.0038630477
f__Rhodobacteraceae 2.924483e-04 0.0038242573
f__Weeksellaceae 2.739210e-04 0.0031292133
f__Eubacteriaceae 1.627561e-04 0.0006691724
f__Sphingobacteriaceae 1.488164e-04 0.0012387571
f__Devosiaceae 1.472568e-04 0.0015006126
f__Pumilibacteraceae 1.262477e-04 0.0007602888
f__WRAU01 9.491039e-05 0.0012411144
f__Peptostreptococcaceae 2.260586e-05 0.0002956099 
family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.2.2 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_rz7gr6doteiz7ad6fyv2
genus mean sd
g__Bacteroides 1.337013e-01 0.0929327655
g__Parabacteroides 9.568074e-02 0.0803866834
g__Phocaeicola 6.854013e-02 0.0793558268
g__Mycoplasmoides 3.015602e-02 0.0749534641
g__Helicobacter_J 2.973811e-02 0.0592188466
g__Odoribacter 2.522456e-02 0.0267412147
g__Roseburia 2.356827e-02 0.0556733585
g__NHYM01 2.223671e-02 0.0792620634
g__Alistipes 2.182523e-02 0.0282922451
g__Coprobacillus 1.990305e-02 0.0873936500
g__Agrobacterium 1.578130e-02 0.0766280394
g__Corynebacterium 1.228390e-02 0.1079320546
g__Akkermansia 1.213003e-02 0.0315247341
g__Fusobacterium_A 1.026752e-02 0.0281722895
g__Kineothrix 8.698172e-03 0.0406735379
g__Proteus 8.556721e-03 0.0677872967
g__Dielma 8.477779e-03 0.0090773537
g__CAG-95 8.014504e-03 0.0203997532
g__JAAYNV01 7.901098e-03 0.0194443338
g__UBA866 7.175409e-03 0.0291873132
g__Desulfovibrio 6.936071e-03 0.0210368286
g__Enterococcus 6.919414e-03 0.0453969341
g__Ureaplasma 6.356785e-03 0.0137373135
g__Lactococcus 6.291980e-03 0.0340225424
g__Lacrimispora 5.993828e-03 0.0097829951
g__Parabacteroides_B 5.924566e-03 0.0099792665
g__CALXRO01 5.698293e-03 0.0306767185
g__Citrobacter 5.650169e-03 0.0332631355
g__NSJ-61 5.495080e-03 0.0197757124
g__Breznakia 5.440147e-03 0.0235692957
g__Clostridium_AQ 5.263895e-03 0.0121113284
g__Salmonella 5.073541e-03 0.0183551010
g__Bilophila 4.925550e-03 0.0088098210
g__Hungatella_A 4.755524e-03 0.0095127458
g__MGBC136627 4.313746e-03 0.0162110012
g__Escherichia 4.139378e-03 0.0264569200
g__UMGS1251 4.111189e-03 0.0072426766
g__Hungatella 4.068243e-03 0.0189714607
g__Clostridium_Q 3.965072e-03 0.0052000165
g__Brevinema 3.962784e-03 0.0122430505
g__Thomasclavelia 3.856935e-03 0.0108480346
g__Scatousia 3.607252e-03 0.0102197530
g__Enterocloster 3.578040e-03 0.0047221634
g__Mailhella 3.569833e-03 0.0101940649
g__Copromonas 3.557270e-03 0.0050035176
g__Ventrimonas 3.472264e-03 0.0070914148
g__Caccovivens 3.304213e-03 0.0121527504
g__Lawsonia 3.250547e-03 0.0116543597
g__Fournierella 3.179598e-03 0.0062039376
g__Limenecus 3.126281e-03 0.0065432982
g__MGBC133411 3.006834e-03 0.0074209199
g__Mucinivorans 2.866176e-03 0.0371005376
g__Sarcina 2.808963e-03 0.0170355922
g__Acetatifactor 2.706113e-03 0.0058179611
g__Eisenbergiella 2.665864e-03 0.0067990748
g__Bacteroides_G 2.651345e-03 0.0344120334
g__CAJLXD01 2.603014e-03 0.0087025932
g__Blautia 2.528781e-03 0.0061107674
g__C-19 2.248901e-03 0.0048466847
g__Butyricimonas 2.191467e-03 0.0049842986
g__Velocimicrobium 2.175479e-03 0.0066412559
g__CAZU01 2.170696e-03 0.0065598332
g__MGBC116941 2.134904e-03 0.0093302936
g__Intestinimonas 2.057548e-03 0.0039214573
g__Negativibacillus 2.044877e-03 0.0054857591
g__Rikenella 1.962168e-03 0.0036910712
g__Phascolarctobacterium 1.921070e-03 0.0050068934
g__RGIG6463 1.768909e-03 0.0039495868
g__JALFVM01 1.721982e-03 0.0038441918
g__Oscillibacter 1.492792e-03 0.0024898538
g__Pseudoflavonifractor 1.485397e-03 0.0027590442
g__Acinetobacter 1.468042e-03 0.0096884882
g__Citrobacter_A 1.376632e-03 0.0060011629
g__Staphylococcus 1.345783e-03 0.0050584305
g__RGIG4733 1.288541e-03 0.0040266511
g__UBA1436 1.199879e-03 0.0064647185
g__Lachnotalea 1.194066e-03 0.0048896855
g__Ruthenibacterium 1.187962e-03 0.0053331515
g__14-2 1.170787e-03 0.0095739735
g__Beduini 1.160220e-03 0.0025025910
g__Scatocola 1.108079e-03 0.0044726561
g__Faecisoma 1.072888e-03 0.0055281036
g__Enterococcus_A 1.071313e-03 0.0098572488
g__Faecimonas 9.764964e-04 0.0082595968
g__RGIG9287 9.526092e-04 0.0092684203
g__CAG-345 9.379861e-04 0.0059852367
g__Blautia_A 9.100332e-04 0.0028927770
g__RGIG1896 8.246024e-04 0.0062408705
g__CAG-269 7.889086e-04 0.0046906408
g__Marseille-P3106 7.845423e-04 0.0017535790
g__WRHT01 6.354364e-04 0.0026829604
g__Eggerthella 6.332937e-04 0.0021152288
g__IOR16 6.324101e-04 0.0020540977
g__Anaerotruncus 6.192165e-04 0.0016862064
g__RUG14156 6.079829e-04 0.0022026464
g__CHH4-2 6.073170e-04 0.0019890673
g__Serratia_A 5.792071e-04 0.0075741155
g__CAG-56 4.857941e-04 0.0016254203
g__Merdimorpha 4.680638e-04 0.0012990360
g__MGBC140009 4.624605e-04 0.0023930648
g__CALURL01 4.581073e-04 0.0016646244
g__Aliarcobacter 4.550604e-04 0.0050029784
g__Scatenecus 4.467347e-04 0.0019650337
g__RGIG8482 4.347613e-04 0.0029582243
g__Enterobacter 4.025794e-04 0.0041076314
g__Klebsiella 4.007019e-04 0.0048624260
g__Caccenecus 3.895102e-04 0.0017702443
g__Alcaligenes 3.656162e-04 0.0047810524
g__Plesiomonas 3.590755e-04 0.0026947991
g__HGM05232 3.575580e-04 0.0009726557
g__JAHHSE01 3.550107e-04 0.0014818046
g__Egerieousia 3.379555e-04 0.0011520742
g__Emergencia 3.373704e-04 0.0017351293
g__Enterococcus_B 3.313107e-04 0.0022138492
g__Stoquefichus 2.990680e-04 0.0020385612
g__Hepatoplasma 2.954146e-04 0.0038630477
g__Paracoccus 2.924483e-04 0.0038242573
g__Moheibacter 2.739210e-04 0.0031292133
g__Scatomorpha 2.610126e-04 0.0010128258
g__UBA7185 2.405856e-04 0.0014474682
g__Eubacterium 1.627561e-04 0.0006691724
g__Sphingobacterium 1.488164e-04 0.0012387571
g__Devosia 1.472568e-04 0.0015006126
g__Anaerosporobacter 1.437105e-04 0.0012673026
g__Caccomorpha 1.366945e-04 0.0010479478
g__UBA2658 1.292159e-04 0.0007163845
g__Protoclostridium 1.262477e-04 0.0007602888
g__Angelakisella 1.253643e-04 0.0009168117
g__Cetobacterium_A 9.201325e-05 0.0008667719
g__Rahnella 6.395024e-05 0.0008362579
g__Peptostreptococcus 2.260586e-05 0.0002956099 
genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.3 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

4.3.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric ) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.4 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

4.5 Permanovas

4.5.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

4.5.1 1. Are the wild populations similar?

4.5.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

4.5.1.2 Number of samples used

[1] 28
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts_func, q=1, dist=dist)

4.5.1.3 Richness

Df SumOfSqs R2 F Pr(>F)
species 1 1.631953 0.207198 6.795072 0.001
Residual 26 6.244347 0.792802 NA NA
Total 27 7.876300 1.000000 NA NA

4.5.1.4 Neutral

Df SumOfSqs R2 F Pr(>F)
species 1 2.018797 0.2566342 8.976049 0.001
Residual 26 5.847641 0.7433658 NA NA
Total 27 7.866438 1.0000000 NA NA

4.5.1.5 Phylogenetic

Df SumOfSqs R2 F Pr(>F)
species 1 0.3786052 0.2108638 6.947419 0.001
Residual 26 1.4168908 0.7891362 NA NA
Total 27 1.7954960 1.0000000 NA NA

4.5.1.6 Functional

Df SumOfSqs R2 F Pr(>F)
species 1 0.0787916 0.1594096 4.930642 0.06
Residual 26 0.4154800 0.8405904 NA NA
Total 27 0.4942716 1.0000000 NA NA

4.5.2 2. Do the antibiotics work?

4.5.2.1 Acclimation vs antibiotics

treat <- meta  %>% #antibiotics samples giving error when plotting
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

4.5.2.2 Number of samples used

[1] 55
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts_func, q=1, dist=dist)
4.5.2.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.042510 0.0911466 6.563244 0.001
species 1 2.300885 0.1026765 7.393486 0.001
individual 26 9.974365 0.4451038 1.232725 0.003
Residual 26 8.091314 0.3610731 NA NA
Total 54 22.409075 1.0000000 NA NA
4.5.2.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.189257 0.1017464 8.979451 0.001
species 1 3.046568 0.1415902 12.495800 0.001
individual 26 9.941985 0.4620568 1.568386 0.001
Residual 26 6.338992 0.2946066 NA NA
Total 54 21.516802 1.0000000 NA NA
4.5.2.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.1642230 0.2000313 18.460841 0.001
species 1 0.8686306 0.0802844 7.409427 0.001
individual 26 4.7385071 0.4379630 1.554596 0.001
Residual 26 3.0480626 0.2817214 NA NA
Total 54 10.8194233 1.0000000 NA NA
4.5.2.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.2832796 0.3550091 33.1502018 0.001
species 1 0.0227598 0.0035387 0.3304425 0.708
individual 26 2.3347730 0.3630154 1.3037622 0.191
Residual 26 1.7907967 0.2784368 NA NA
Total 54 6.4316092 1.0000000 NA NA 
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))


beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

4.5.3 3. Do the FMT work?

4.5.3.1 Comparison between FMT2 vs Post-FMT2

transplant3 <- meta  %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))

transplant3.counts_func <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts_func)),sort(as.character(rownames(transplant3))))

transplant3_nmds <- sample_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

4.5.3.2 Number of samples used

[1] 54
beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts_func, q=1, dist=dist)
4.5.3.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.9748558 0.0733017 5.103443 0.001
time_point 1 0.7043146 0.0529590 3.687140 0.001
type 1 1.4084037 0.1059011 7.373099 0.001
individual 25 7.1553504 0.5380273 1.498352 0.001
Residual 16 3.0563076 0.2298108 NA NA
Total 44 13.2992322 1.0000000 NA NA
4.5.3.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 1.1035252 0.0842380 6.734701 0.001
time_point 1 0.7327432 0.0559342 4.471856 0.001
type 1 1.8565350 0.1417193 11.330241 0.001
individual 25 6.7855728 0.5179794 1.656466 0.001
Residual 16 2.6217058 0.2001290 NA NA
Total 44 13.1000821 1.0000000 NA NA
4.5.3.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1387836 0.0608025 2.904615 0.126
time_point 1 0.1416869 0.0620745 2.965380 0.002
type 1 0.1763572 0.0772639 3.690999 0.002
individual 25 1.0612162 0.4649298 0.888412 0.393
Residual 16 0.7644858 0.3349292 NA NA
Total 44 2.2825297 1.0000000 NA NA
4.5.3.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 -0.0003437 -0.0001753 -0.0062959 0.903
time_point 1 0.0253408 0.0129285 0.4642107 0.713
type 1 0.1178144 0.0601071 2.1582044 0.137
individual 25 0.9438371 0.4815314 0.6915940 0.747
Residual 16 0.8734253 0.4456084 NA NA
Total 44 1.9600740 1.0000000 NA NA 
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

4.5.4 4. Are there differences between the control and the treatment group?

4.5.4.1 Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2.counts_func <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts_func)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

4.5.4.2 Number of samples used

[1] 28
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts_func, q=1, dist=dist)
4.5.4.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.8966107 0.1132233 3.515588 0.001
type 1 0.6463814 0.0816245 2.534445 0.001
Residual 25 6.3759661 0.8051521 NA NA
Total 27 7.9189582 1.0000000 NA NA
4.5.4.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 0.9398968 0.1224370 4.112304 0.001
type 1 1.0227481 0.1332297 4.474800 0.001
Residual 25 5.7139316 0.7443333 NA NA
Total 27 7.6765766 1.0000000 NA NA
4.5.4.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1200820 0.1454923 4.647187 0.001
type 1 0.0592745 0.0718175 2.293931 0.031
Residual 25 0.6459931 0.7826902 NA NA
Total 27 0.8253496 1.0000000 NA NA
4.5.4.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 0.0116738 0.0168395 0.4229082 0.530
type 1 -0.0085273 -0.0123007 -0.3089208 0.915
Residual 25 0.6900904 0.9954612 NA NA
Total 27 0.6932368 1.0000000 NA NA 
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
grid.arrange(arrangeGrob(p0,p1,p2,p3))